source("~/jules_ppe_gmd/docs/JULES-ES-1p0-common-packages.R")
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.7-0 (2021-06-25) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: viridis
## Loading required package: viridisLite
## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code
## Loading required package: lhs
## Loading required package: parallel
## Loading required package: viztools
source("~/jules_ppe_gmd/docs/JULES-ES-1p0-common-functions.R")
## ----------------------------------------------------------------------
## Data locations and constants
## ----------------------------------------------------------------------
#ensloc <- '/project/carbon_ppe/JULES-ES-1p0_PPE/'
ensloc_wave00 <- '/data/users/hadaw/JULES_ES_PPE/u-au932/'
ensloc_wave01 <- '/data/users/hadaw/JULES_ES_PPE/u-ck006/'
# Some pallete options
yg = brewer.pal(9, "YlGn")
ryb = brewer.pal(11, "RdYlBu")
byr = rev(ryb)
rb = brewer.pal(11, "RdBu")
br = rev(rb)
blues = brewer.pal(9, 'Blues')
cbPal <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
wave00col <- 'skyblue2'
wave01col <- 'tomato2'
zissou5 <- wes_palette('Zissou1', 5, type = c('discrete', 'continuous'))
zblue <- makeTransparent(as.character(zissou5)[1], 150)
zred <- makeTransparent(as.character(zissou5)[5], 150)
ysec = 60*60*24*365
years <- 1850:2013
global_mean_modern_value_wave00_file <- "~/jules_ppe_gmd/docs/data/global_mean_modern_value_wave00_2022-13-09.rdata"
load(global_mean_modern_value_wave00_file)
modern_value_stan_file <- "~/jules_ppe_gmd/docs/data/modern_value_stan_2022-09-13.rdata"
load(modern_value_stan_file)
Y <- datmat_wave00
Y_nlevel0_ix <- which(is.na(datmat_wave00[,'year']))
# Load up the data
lhs_i = read.table('~/jules_ppe_gmd/docs/data/lhs_u-ao732.txt', header = TRUE)
lhs_ii = read.table('~/jules_ppe_gmd/docs/data/lhs_u-ao732a.txt', header = TRUE)
toplevel_ix = 1:499
# The raw input data is a latin hypercube
lhs = rbind(lhs_i, lhs_ii)[toplevel_ix, ]
lhs_level0 <- lhs[-Y_nlevel0_ix,]
X = normalize(lhs)
colnames(X) = colnames(lhs)
X_level0 <- X[-Y_nlevel0_ix,]
X_nlevel0 <- X[Y_nlevel0_ix,]
d = ncol(X)
# lower and higher bound on the normalised matrix for visualisation
rx = rbind(rep(0,32), rep(1,32))
low_npp_ix <- which(Y[,'npp_nlim_lnd_sum'] < 1e5)
# code from https://stackoverflow.com/questions/28182872/how-to-use-different-sets-of-data-in-lower-and-upper-panel-of-pairs-function-in
#X <- matrix(runif(300), ncol=3)
#Y <- matrix(c(sort(runif(100, 0, 10)),
# sort(runif(100, 0, 10)),
# sort(runif(100, 0, 10))), ncol=3)
#pdf(file = 'figs/fig02.pdf', width = 12, height = 10)
#pdf(file = 'figs/run-failure-pairs.pdf', width = 12, height = 10)
x1 <- X[low_npp_ix, ]
x2 <- X_nlevel0
XY <- rbind(x1, x2)
pairs(XY,
lower.panel=function(x, y, ...) {
Xx <- x[seq_len(nrow(x1))] # corresponds to X subset
Xy <- y[seq_len(nrow(x1))] # corresponds to X subset
#usr <- par("usr"); on.exit(par(usr))
#par(usr = c(range(x1[, -ncol(x1)]), range(x1[, -1]))) # set up limits
points(Xx, Xy, col = zblue, pch = 19, cex = 0.8)
# if(par('mfg')[2] == 1) axis(2) # if left plot, add left axis
#if(par('mfg')[1] == ncol(x1)) axis(1) # if bottom plot add bottom axis
},
upper.panel=function(x, y, ...) {
Yx <- x[(nrow(x1) + 1):length(x)] # Y subset
Yy <- y[(nrow(x1) + 1):length(y)] # Y subset
#cntr <- outer(Yx, Yx, FUN='*') # arbitrary function for contour
# usr <- par("usr"); on.exit(par(usr))
#par(usr = c(range(x2[, -1]), range(x2[, -ncol(x2)]))) # set up limits
points(Yx, Yy, col = zred, pch = 19, cex = 0.8)
#contour(Yx, Yy, cntr, add=TRUE)
#if(par('mfg')[2] == ncol(x2)) axis(4) # if right plot, add right axis
#if(par('mfg')[1] == 1) axis(3) # if top plot, add top axis
},
#tick=FALSE, # suppress the default tick marks
#line=1,
gap = 0,
xlim = c(0,1), ylim = c(0,1),
labels = 1:d,
oma = c(2, 18, 2, 2)) # move the default tick labels off the plot
reset()
legend('left', legend = paste(1:d, colnames(lhs)), cex = 1.1, bty = 'n')
legend('topleft', pch = 19, col = c( zred, zblue), legend = c('failed', 'zero carbon cycle'), bty = 'n', inset = 0.02, cex = 1.1 )

#dev.off()
Generate the data to save
Y_char <- rep('RAN', nrow(X))
Y_char[Y_nlevel0_ix] <- 'CRASHED'
Y_char[low_npp_ix] <- 'LOWCARBON'
Y_factor <- as.factor(Y_char)
wave0_summary_raw <- cbind(lhs, Y, Y_factor)
wave0_summary_df <- cbind(X, Y, Y_factor)
save(lhs, X, Y, Y_char, Y_factor, wave0_summary_raw, wave0_summary_df, file = 'data/wave0_summary.RData')